#include "wrapper_sstructmg.h"

struct OrigBlkInfo {
    int pid = -1;
    idx_t coords[2] = {-1, -1};
    idx_t begs[2] = {-1, -1};
    idx_t ends[2] = {-1, -1};
};

struct MergBlkInfo {
    idx_t bid = -1;
    int nprocs = 0;
    idx_t begs[2] = {INT_MAX, INT_MAX};
    idx_t ends[2] = {INT_MIN, INT_MIN};
    bool operator < (const MergBlkInfo & b) const {
        if (begs[1] == b.begs[1]) {
            if (nprocs == b.nprocs) return begs[0] < b.begs[0];
            else return nprocs < b.nprocs;
        }
        else return begs[1] < b.begs[1];
    }
};

const idx_t num_diag = 5;

int main(int argc, char* argv[])
{
    int my_pid, num_procs;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_pid);
    MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
    if (my_pid == 0) printf("\033[1;35mMax threads %d\033[0m \n", omp_get_max_threads());

    idx_t arg_cnt = 1;
    const idx_t glb_dims [2] = { atoi(argv[arg_cnt++]), atoi(argv[arg_cnt++]) };
    const idx_t asg_dims [2] = { atoi(argv[arg_cnt++]), atoi(argv[arg_cnt++]) };// assigned block dimensions
    const std::string pathname(argv[arg_cnt++]);
    const std::string   its_name = std::string(argv[arg_cnt++]),
                        prc_name = std::string(argv[arg_cnt++]);
    std::string config_mg_file;
	if (prc_name == "GMG") config_mg_file = std::string(argv[arg_cnt++]);

    if (my_pid == 0) printf(" Glb %d %d Blk %d %d\n", glb_dims[0], glb_dims[1], asg_dims[0], asg_dims[1]);
    MPI_Barrier(MPI_COMM_WORLD);

    std::vector<OrigBlkInfo> orig_blk_infos;

    const std::string distinfo_name = "distinfo." + std::to_string(glb_dims[1]) + "x" + std::to_string(glb_dims[0])
                                        + ".tile" + std::to_string(asg_dims[1]) + "x" + std::to_string(asg_dims[0]);
    {
        FILE * fp = fopen((pathname + "/" + distinfo_name).c_str(), "r"); assert(fp);
        fscanf(fp, "i1 i2 j1 j2 i j Index_withLand Index_rmvland");
        const idx_t num_lines = (glb_dims[0] / asg_dims[0]) * (glb_dims[1] / asg_dims[1]);
        for (idx_t b = 0; b < num_lines; b++) {
            idx_t i1, i2, j1, j2, pi, pj, index_withLand, index_rmvLand;
            fscanf(fp, "%d,%d,%d,%d,%d,%d,%d,%d", & i1, & i2, & j1, & j2, & pi, & pj, & index_withLand, & index_rmvLand);
            if (index_rmvLand == -9999) continue;

            OrigBlkInfo blk;
            blk.pid = orig_blk_infos.size(); assert(blk.pid == index_rmvLand - 1);
            blk.coords[1] = pi - 1;
            blk.coords[0] = pj - 1;

            blk.begs[1] = i1 - 1;
            blk.ends[1] = i2; assert(blk.ends[0] <= glb_dims[0]);

            blk.begs[0] = j1 - 1;
            blk.ends[0] = j2; assert(blk.ends[1] <= glb_dims[1]);

            assert(blk.ends[0] - blk.begs[0] == asg_dims[0]);
            assert(blk.ends[1] - blk.begs[1] == asg_dims[1]);
            orig_blk_infos.emplace_back(blk);
        }
        fclose(fp);
    }
    const idx_t orig_nblks = orig_blk_infos.size(); assert(orig_nblks == num_procs);

    // 重组合并一些原来的块
    std::vector<idx_t> merged_blk_ptr;
    for (idx_t j = 0; j < orig_nblks; j++) {
        if (j == 0 ||   orig_blk_infos[j].begs[0] != orig_blk_infos[j-1].begs[0] ||
                        orig_blk_infos[j].begs[1] != orig_blk_infos[j-1].ends[1] )
        {
            merged_blk_ptr.push_back(j);
        } else {
            assert(orig_blk_infos[j].ends[0] == orig_blk_infos[j-1].ends[0]);
        }
    }
    merged_blk_ptr.push_back(orig_nblks);

    const idx_t merged_nblks = merged_blk_ptr.size() - 1;
    std::vector<idx_t> map_orig2merged(orig_nblks, -1);
    std::vector<MergBlkInfo> merged_blk_infos;
    for (idx_t mb = 0; mb < merged_nblks; mb++) {
        MergBlkInfo mbInfo;
        mbInfo.bid = mb;
        mbInfo.nprocs = merged_blk_ptr[mb + 1] - merged_blk_ptr[mb];
        for (idx_t p = merged_blk_ptr[mb]; p < merged_blk_ptr[mb + 1]; p++) {
            map_orig2merged[p] = mb;

            mbInfo.begs[0] = std::min(mbInfo.begs[0], orig_blk_infos[p].begs[0]);
            mbInfo.begs[1] = std::min(mbInfo.begs[1], orig_blk_infos[p].begs[1]);

            mbInfo.ends[0] = std::max(mbInfo.ends[0], orig_blk_infos[p].ends[0]);
            mbInfo.ends[1] = std::max(mbInfo.ends[1], orig_blk_infos[p].ends[1]);
        }
        merged_blk_infos.push_back(mbInfo);
    }
#ifdef DEBUG_PRINT
    if (my_pid == 0) {
        for (idx_t mb = 0; mb < merged_nblks; mb++) {
            const MergBlkInfo & mbInfo = merged_blk_infos[mb];
            printf("Merged Blk %d : (%d,%d) ~ (%d,%d) Proc ", mb, mbInfo.begs[0], mbInfo.begs[1], mbInfo.ends[0], mbInfo.ends[1]);
            for (idx_t p = merged_blk_ptr[mb]; p < merged_blk_ptr[mb + 1]; p++) {
                printf("%d ", p);
            }
            printf("\n");
        }
    }
    MPI_Barrier(MPI_COMM_WORLD);
#endif
#if 0
    {// stage 2: dimension order merging
        std::vector<MergBlkInfo> tmp_arr = merged_blk_infos;
        std::sort(tmp_arr.begin(), tmp_arr.end());

        // if (my_pid == 0) {
        //     for (idx_t mb = 0; mb < merged_nblks; mb++) {
        //         const MergBlkInfo & mbInfo = tmp_arr[mb];
        //         printf("Sorted Blk %d : (%d,%d) ~ (%d,%d) Bid %d #proc %d",
        //             mb, mbInfo.begs[0], mbInfo.begs[1], mbInfo.ends[0], mbInfo.ends[1], mbInfo.bid, mbInfo.nprocs);
        //         printf("\n");
        //     }
        // }
        // MPI_Barrier(MPI_COMM_WORLD);

        std::vector<idx_t> dom_cnt;
        std::vector<idx_t> dom_ids;

        idx_t out_ptr0 = 0, out_ptr1 = 0;
        while (out_ptr0 < merged_nblks) {
            while (out_ptr1 < merged_nblks  && tmp_arr[out_ptr1].nprocs == tmp_arr[out_ptr0].nprocs
                                            && tmp_arr[out_ptr1].begs[1]== tmp_arr[out_ptr0].begs[1])
            {// 同长度且起始位置相同的才有可能合并
                out_ptr1 ++;
            }

            // if (my_pid == 0) printf("Out %d %d\n", out_ptr0, out_ptr1);

            // 在 [out_ptr0, out_ptr1) 内寻找可以合并的
            // 继续使用双指针搜索
            idx_t in_ptr0 = out_ptr0, in_ptr1 = in_ptr0;
            while (in_ptr0 < out_ptr1) {
                idx_t last_e0 = tmp_arr[in_ptr0].ends[0];
                while (in_ptr1 < out_ptr1 && (in_ptr0 == in_ptr1 || last_e0 == tmp_arr[in_ptr1].begs[0]) ) {
                    last_e0 = tmp_arr[in_ptr1].ends[0];// 更新
                    in_ptr1 ++;
                }

                // if (my_pid == 0) printf("  In %d %d\n", in_ptr0, in_ptr1);

                // 跳出循环时
                if (in_ptr1 > in_ptr0) {
                    dom_cnt.push_back(in_ptr1 - in_ptr0);
                    for (idx_t p = in_ptr0; p < in_ptr1; p++) dom_ids.push_back(tmp_arr[p].bid);
                }
                in_ptr0 = in_ptr1;// 移动指针
            }


            out_ptr0 = out_ptr1;// 移动指针
        }

        if (my_pid == 0) {
            idx_t off = 0;
            for (idx_t domb = 0; domb < (idx_t) dom_cnt.size(); domb++) {
                printf("Dom_Blk %d cnt %d: Merged_Blk ", domb, dom_cnt[domb]);
                for (idx_t p = off; p < off + dom_cnt[domb]; p++) printf("%d ", dom_ids[p]);
                printf("\n");
                off += dom_cnt[domb];
            } assert(off == merged_nblks);
        }
        MPI_Barrier(MPI_COMM_WORLD);
    }
#endif
    const OrigBlkInfo & my_orig_blk = orig_blk_infos[my_pid];
    const idx_t my_part = map_orig2merged[my_pid]; assert(my_part < merged_nblks);// 本进程属于哪个块
    const MergBlkInfo & my_merg_blk = merged_blk_infos[my_part];

    const idx_t * glb_begs = my_orig_blk.begs, * glb_ends = my_orig_blk.ends;// 全局的起始和结束位置
    idx_t box_begs[3], box_ends[3];// 本进程所负责区域在融合后的块内的位置

    box_begs[0] = 0;
    box_ends[0] = 1;

    box_begs[1] = glb_begs[0] - my_merg_blk.begs[0];
    box_ends[1] = box_begs[1] + (glb_ends[0] - glb_begs[0]);

    box_begs[2] = glb_begs[1] - my_merg_blk.begs[1];
    box_ends[2] = box_begs[2] + (glb_ends[1] - glb_begs[1]);

    const idx_t blk_exts[3] = {1,   my_merg_blk.ends[0] - my_merg_blk.begs[0],
                                    my_merg_blk.ends[1] - my_merg_blk.begs[1] };

    const idx_t my_nelems = (glb_ends[0] - glb_begs[0]) * (glb_ends[1] - glb_begs[1]);
    assert(my_nelems == (box_ends[0] - box_begs[0]) * (box_ends[1] - box_begs[1]) * (box_ends[2] - box_begs[2]));

    // 先读入本进程所负责的区域的数据
    data_t * A_buf = new data_t [my_nelems * num_diag];
    data_t * rhs_buf = new data_t [my_nelems];// 读入数据时从外到内为 [3360][3072]
    data_t * x0_buf  = new data_t [my_nelems];// 读入数据时从外到内为 [3360][3072]

    {// 读入向量
        MPI_File fh = MPI_FILE_NULL;// 文件句柄
        MPI_Datatype etype = MPI_DOUBLE;
        MPI_Datatype read_type = MPI_DATATYPE_NULL;// 写出类型
        idx_t size[2], subsize[2], start[2];
        size   [0] = glb_dims[0]; size   [1] = glb_dims[1];
        subsize[0] = asg_dims[0]; subsize[1] = asg_dims[1];
        start  [0] = glb_begs[0]; start  [1] = glb_begs[1];

        MPI_Type_create_subarray(2, size, subsize, start, MPI_ORDER_C, etype, &read_type);
        MPI_Type_commit(&read_type);
        MPI_Offset displacement = 0;
        displacement *= sizeof(double);// 位移要以字节为单位！
        MPI_Status status;

        int ret;
        ret = MPI_File_open(MPI_COMM_WORLD, (pathname + "/array_b").c_str(), MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);
        if (ret != MPI_SUCCESS) {
            printf("Could not open b file: ret %d\n", ret);
        }
        MPI_File_set_view(fh, displacement, etype, read_type, "native", MPI_INFO_NULL);
        MPI_File_read_all(fh, rhs_buf, my_nelems, etype, &status);
        MPI_File_close(&fh);

        ret = MPI_File_open(MPI_COMM_WORLD, (pathname + "/array_x0").c_str(), MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);
        if (ret != MPI_SUCCESS) {
            printf("Could not open x file: ret %d\n", ret);
        }
        MPI_File_set_view(fh, displacement, etype, read_type, "native", MPI_INFO_NULL);
        MPI_File_read_all(fh, x0_buf, my_nelems, etype, &status);
        MPI_File_close(&fh);

        MPI_Type_free(&read_type);
    }
    {// 读入矩阵
        MPI_File fh = MPI_FILE_NULL;// 文件句柄
        MPI_Datatype etype = MPI_DOUBLE;
        MPI_Datatype read_type = MPI_DATATYPE_NULL;// 写出类型
        idx_t size[3], subsize[3], start[3];
        // 注意这里要反过来，因为存储的数据从外到内是[3072][3360]的
        size   [0] = glb_dims[0]; size   [1] = glb_dims[1]; size   [2] = num_diag;
        subsize[0] = asg_dims[0]; subsize[1] = asg_dims[1]; subsize[2] = num_diag;
        start  [0] = glb_begs[0]; start  [1] = glb_begs[1]; start  [2] = 0;

        MPI_Type_create_subarray(3, size, subsize, start, MPI_ORDER_C, etype, &read_type);
        MPI_Type_commit(&read_type);
        MPI_Offset displacement = 0;
        displacement *= sizeof(double);// 位移要以字节为单位！
        MPI_Status status;

        int ret;
        ret = MPI_File_open(MPI_COMM_WORLD, (pathname + "/mat.AOS.bin").c_str(), MPI_MODE_RDONLY, MPI_INFO_NULL, &fh);
        if (ret != MPI_SUCCESS) {
            printf("Could not open file: ret %d\n", ret);
        }
        MPI_File_set_view(fh, displacement, etype, read_type, "native", MPI_INFO_NULL);
        MPI_File_read_all(fh, A_buf, my_nelems * num_diag, etype, &status);
        MPI_File_close(&fh);
        MPI_Type_free(&read_type);
    }
    {// 手动检验
        double loc_b_dot = 0.0, loc_x_dot = 0.0, loc_A_dot = 0.0;
        for (idx_t i = 0; i < my_nelems; i++) {
            loc_b_dot += rhs_buf[i] * rhs_buf[i];
            loc_x_dot +=  x0_buf[i] *  x0_buf[i];
            if (A_buf[i * num_diag + 2] == -1.0) continue;
            // for (idx_t d = 0; d < num_diag; d++)
                loc_A_dot += A_buf[i * num_diag + 2] * A_buf[i * num_diag + 2];
        }
        // printf("Proc %d : %.15e %.15e\n", my_pid, loc_b_dot, loc_x_dot);
        double ck_b_dot = 0.0, ck_x_dot = 0.0, ck_A_dot = 0.0;
        MPI_Allreduce(& loc_b_dot, & ck_b_dot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        MPI_Allreduce(& loc_x_dot, & ck_x_dot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        MPI_Allreduce(& loc_A_dot, & ck_A_dot, 1, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
        if (my_pid == 0) {
            printf("My Calc:\n");
            printf("(  b,   b) = %.20e\n",  (double)ck_b_dot);
            printf("(  x,   x) = %.20e\n",  (double)ck_x_dot);
            printf("(  A,   A) = %.20e\n",  (double)ck_A_dot);
        }
    }
    
    std::vector< Entry<idx_t, data_t> > all_nnz;
    {
        // 确定自己的四个方向的邻居块号
        idx_t i1_n_pid = -1, i1_p_pid = -1, i2_n_pid = -1, i2_p_pid = -1;

        for (size_t k = 0; k < orig_blk_infos.size(); k++) {
            const OrigBlkInfo & ngb_blk = orig_blk_infos[k];
            if (ngb_blk.begs[0] == glb_begs[0] && ngb_blk.ends[0] == glb_ends[0]) {// i1维相同
                if (ngb_blk.begs[1] == glb_ends[1]) {
                    assert(i2_p_pid == -1);
                    i2_p_pid = ngb_blk.pid;
                }
                if (ngb_blk.ends[1] == glb_begs[1]) {
                    assert(i2_n_pid == -1);
                    i2_n_pid = ngb_blk.pid;
                }
            }
            if (ngb_blk.begs[1] == glb_begs[1] && ngb_blk.ends[1] == glb_ends[1]) {// i2维相同
                if (ngb_blk.begs[0] == glb_ends[0]) {
                    assert(i1_p_pid == -1);
                    i1_p_pid = ngb_blk.pid;
                }
                if (ngb_blk.ends[0] == glb_begs[0]) {
                    assert(i1_n_pid == -1);
                    i1_n_pid = ngb_blk.pid;
                }
            }
        }
        // if (i1_n_pid != -1 && map_orig2merged[i1_n_pid] == my_part) i1_n_pid = -1;
        // if (i1_p_pid != -1 && map_orig2merged[i1_p_pid] == my_part) i1_p_pid = -1;
        // if (i2_n_pid != -1 && map_orig2merged[i2_n_pid] == my_part) i2_n_pid = -1;
        // if (i2_p_pid != -1 && map_orig2merged[i2_p_pid] == my_part) i2_p_pid = -1;
#ifdef DEBUG_PRINT
        for (int p = 0; p < num_procs; p++) {
            if (p == my_pid) {
                printf("Proc %3d i1 %3d %3d i2 %3d %3d Box (%5d,%5d,%5d)~(%5d,%5d,%5d) Blk %3d (%5d,%5d,%5d)\n",
                    my_pid, i1_n_pid, i1_p_pid, i2_n_pid, i2_p_pid, 
                    box_begs[0], box_begs[1], box_begs[2], box_ends[0], box_ends[1], box_ends[2], my_part, blk_exts[0], blk_exts[1], blk_exts[2]);
                fflush(stdout);
            }
            MPI_Barrier(MPI_COMM_WORLD);
        }
#endif
        // 提取块边界非零元
        #pragma omp parallel for collapse(2) schedule(static)
        for (idx_t i0 = box_begs[0]; i0 < box_ends[0]; i0++)
        for (idx_t i1 = box_begs[1]; i1 < box_ends[1]; i1++)
        for (idx_t i2 = box_begs[2]; i2 < box_ends[2]; i2++) {
            StructDescriptor<idx_t> row, col;
            data_t* dst_ptr = A_buf + (((i0 - box_begs[0] ) * (box_ends[1] - box_begs[1])
                                    +    i1 - box_begs[1] ) * (box_ends[2] - box_begs[2])
                                    +    i2 - box_begs[2] ) * num_diag;
            assert(i0 == 0);
            row.id = my_part;
            row.index[0] = i0; row.index[1] = i1; row.index[2] = i2;
            col.index[0] = i0;

            const idx_t glb_pos[3] = { i0,  i1 + my_merg_blk.begs[0],
                                            i2 + my_merg_blk.begs[1] };// 该点在全局的位置

            if (i1 == 0               && i1_n_pid != -1) {
                col.id = map_orig2merged[i1_n_pid]; 
                // if (col.id == row.id) {
                //     printf("Error !! Proc %d row.id %d i1_n_pid %d col.id %d\n", my_pid, row.id, i1_n_pid, col.id);
                // }
                assert(col.id != row.id);
                const MergBlkInfo & ngb_blk = merged_blk_infos[col.id];
                col.index[1] = glb_pos[1] - 1 - ngb_blk.begs[0];// col.index[1] = blk_exts[1] - 1;
                col.index[2] = glb_pos[2]     - ngb_blk.begs[1];

                // if (my_pid == 0) printf("i1n Row %d (%d,%d,%d) Col %d (%d,%d,%d)\n", row.id, row.index[0], row.index[1], row.index[2],
                //     col.id, col.index[0], col.index[1], col.index[2]);

                all_nnz.emplace_back(row, col, dst_ptr[0]);
                dst_ptr[0] = 0.0;
            }
            if (i1 == blk_exts[1] - 1 && i1_p_pid != -1) {
                col.id = map_orig2merged[i1_p_pid];
                // if (col.id == row.id) {
                //     printf("Error !! Proc %d row.id %d i1_p_pid %d col.id %d\n", my_pid, row.id, i1_p_pid, col.id);
                // }
                assert(col.id != row.id);
                const MergBlkInfo & ngb_blk = merged_blk_infos[col.id];
                col.index[1] = glb_pos[1] + 1 - ngb_blk.begs[0];// col.index[1] = 0;
                col.index[2] = glb_pos[2]     - ngb_blk.begs[1];

                // if (my_pid == 0) printf("i1p Row %d (%d,%d,%d) Col %d (%d,%d,%d)\n", row.id, row.index[0], row.index[1], row.index[2],
                //     col.id, col.index[0], col.index[1], col.index[2]);

                all_nnz.emplace_back(row, col, dst_ptr[4]);
                dst_ptr[4] = 0.0;
            }
            if (i2 == 0               && i2_n_pid != -1) { assert(false);
                col.id = map_orig2merged[i2_n_pid];// col.id = i2_n_bid;
                const MergBlkInfo & ngb_blk = merged_blk_infos[col.id];
                col.index[1] = glb_pos[1]     - ngb_blk.begs[0];
                col.index[2] = glb_pos[2] - 1 - ngb_blk.begs[1];// col.index[2] = blk_exts[2] - 1;
                
                all_nnz.emplace_back(row, col, dst_ptr[1]);
                dst_ptr[1] = 0.0;
            }
            if (i2 == blk_exts[2] - 1 && i2_p_pid != -1) { assert(false);
                col.id = map_orig2merged[i2_p_pid];// col.id = i2_p_bid;
                const MergBlkInfo & ngb_blk = merged_blk_infos[col.id];
                col.index[1] = glb_pos[1]     - ngb_blk.begs[0];
                col.index[2] = glb_pos[2] + 1 - ngb_blk.begs[1];// col.index[2] = 0;
                
                all_nnz.emplace_back(row, col, dst_ptr[3]);
                dst_ptr[3] = 0.0;
            }
        }
    }

    SStructMG_Init();
    {
        const idx_t glb_nparts = merged_nblks;
        SStructGrid<idx_t> ssgrid(MPI_COMM_WORLD, 2, glb_nparts);
        if (my_pid == 0) printf("Total Blks %d\n", glb_nparts);
        const idx_t halos[3] = { 0, 1, 1 };
        
        {// 添加本进程所拥有的块和盒子
            const bool periods[3] = { false, false, false };
            ssgrid.addBlock(my_part, blk_exts, periods);
            ssgrid.addBox(my_part, box_begs, box_ends);
        }
        ssgrid.assemble();

        const idx_t my_nblks = ssgrid.my_blk_list.size(); assert(my_nblks == 1);
        idx_t shalos [my_nblks * 3];
        for (idx_t k = 0; k < my_nblks; k++) {
            shalos[k * 3 + 0] = halos[0];
            shalos[k * 3 + 1] = halos[1];
            shalos[k * 3 + 2] = halos[2];
        }

        // 根据半结构网格建立向量
        par_SstructVector<idx_t, data_t> b(ssgrid, shalos, true);
        par_SstructVector<idx_t, data_t> x(b), y(b), x0(b);
        {// 填入向量数据
            b.set_box_values(my_part, box_begs, box_ends, rhs_buf);
            x0.set_val(0.0, true);
            y.set_val(0.0, true);
        }
        
        // 根据半结构网格建立矩阵
        __int128_t masks[my_nblks]; for (idx_t k = 0; k < my_nblks; k++) masks[k] = stencil_mask_2d5i1i2;
        par_SstructMatrix<idx_t, data_t, data_t> A(ssgrid, shalos, masks);
        {// 填入矩阵数据
            A.init_offd();// 初始化非对角部分
            A.set_box_values(my_part, box_begs, box_ends, A_buf);
            // 然后处理非对角快，也即非结构部分
            A.add_entries(all_nnz.size(), all_nnz.data());
        }
        A.assemble();

        const int test_cnt = 10;
        std::vector<TEST_RECORD> records;
        for (int test = 0; test < test_cnt; test++) {
            vec_copy(x0, x);

            double ck = vec_dot<idx_t, data_t, double>(b, b);
            if (my_pid == 0) printf("(  b,   b) = %.15e\n", ck);

            ck = vec_dot<idx_t, data_t, double>(x, x);
            if (my_pid == 0) printf("(  x,   x) = %.15e\n", ck);

            A.Mult(b, y, false);
            ck = vec_dot<idx_t, data_t, double>(y, y);
            if (my_pid == 0) printf("(A*b, A*b) = %.15e\n", ck);

            A.Mult(x, y, false);
            ck = vec_dot<idx_t, data_t, double>(y, y);
            if (my_pid == 0) printf("(A*x, A*x) = %.15e\n", ck);

            
            TEST_CONFIG config;
            config.my_nblks = my_nblks;
            config.config_mg_file = config_mg_file;
            config.ssgrid_ptr = & ssgrid;
            config.rtol = 1.0e-7;
            TEST_RECORD rec;

            buildup_solver(its_name, prc_name, config);
            setup_and_solve(A, b, x, rec);

            data_t true_r_norm, b_norm;
            check_residual(A, x, b, y, true_r_norm, b_norm);

            if (my_pid == 0) {
                printf("\033[1;35mtrue ||r|| = %20.16e ||r||/||b||= %20.16e\033[0m\n", 
                    true_r_norm, true_r_norm / b_norm);
                printf("Proc %d Setup, Solve costs %.6f %.6f s\n",
                    my_pid, rec.setup, rec.solve);
            }

            stat_part_times(rec);
            records.push_back(rec);

            destroy_solver();
            
        }// test loop

        // 输出最佳时间
        if (my_pid == 0) {
            stat_multi_runs(records);
        }
    }// all Semi-** object must be destroyed before MPI_Finalize()
    SStructMG_Finalize();
    delete [] A_buf;
    delete [] rhs_buf;
    delete [] x0_buf;
    
    MPI_Finalize();
    
    return 0;
}